#####
# Replication for: "Can Political Speech Foster Tolerance of Immigrants?" by Schleiter, Tavits, and Ward.
# Figure S.8
#####

library(here)
library(data.table)
library(lfe)
library(ggplot2)
library(ggstance)
library(cowplot)

#####-----#####-----#####-----
# NOTE:
# Changes to either R or lfe (or lfe's dependencies) since the original fitting in July 2019
# mean that the original code no longer exactly replicates the clustered standard errors
# for the behavioral model. The differences are small and do not affect any substantive
# interpretations.
# It is hard to say what exactly is causing this, although it may be the dof used. It is important
# to note that the coefficient estimates and number of observations still exactly match.
#####-----#####-----#####-----

# get functions (need extractPlotInfo)
source(here("functions.R"))

if(!exists("pilot1")){
  pilot1 <- fread(here("data", "pilot1.csv"))
}

if(!exists("pilot1_behav")){
  pilot1_behav <- fread(here("data", "pilot1_behav.csv"))
}

if(!exists("pilot2")){
  pilot2 <- fread(here("data", "pilot2.csv"))
}


# Fitting -----------------------------------------------------------------

# Pilot 1, attitudinal models

p1_neighbors <- lm(
  imm_neighbors ~ CH + CS + NM + gender + edu + race + news + party_id + age + ageSq + region,
  data = pilot1)

p1_increase <- lm(
  imm_increase ~ CH + CS + NM + gender + edu + race + news + party_id + age + ageSq + region,
  data = pilot1)

p1_safety <- lm(
  imm_safety ~ CH + CS + NM + gender + edu + race + news + party_id + age + ageSq + region,
  data = pilot1)

p1_culture <- lm(
  imm_culture ~ CH + CS + NM + gender + edu + race + news + party_id + age + ageSq + region,
  data = pilot1)

p1_services <- lm(
  imm_services ~ CH + CS + NM + gender + edu + race + news + party_id + age + ageSq + region,
  data = pilot1)

p1_pca_all <- lm(
  immPCA_all ~ CH + CS + NM + gender + edu + race + news + party_id + age + ageSq + region,
  data = pilot1)

# Pilot 1, behavioral model

game_mod <- felm(
  tokens ~ female_partner + foreign_partner + foreign_partner:CH + foreign_partner:CS + foreign_partner:NM | ResponseId + game_NR | 0 | ResponseId,
  data = pilot1_behav)

# Pilot 2, attitudinal models

p2_neighbors <- lm(
  imm_neighbors ~ CH + CS + NM + gender + race + edu + news + age + ageSq + party_id + region,
  data = pilot2)

p2_culture <- lm(
  imm_culture ~ CH + CS + NM + gender + race + edu + news + age + ageSq + party_id + region,
  data = pilot2)

p2_pca <- lm(
  immPCA ~ CH + CS + NM + gender + race + edu + news + age + ageSq + party_id + region,
  data = pilot2)


# Plot preparation --------------------------------------------------------

plot_mat_p1 <- data.frame(est = rep(0, 18), lower95 = 0, upper95 = 0, lower90 = 0, upper90 = 0, x = 0, D = NA, stringsAsFactors = F)
plot_mat_p1[ 4:6, ] <- data.frame(extractPlotInfo(p1_neighbors), 5, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p1[ 7:9, ] <- data.frame(extractPlotInfo(p1_increase), 4, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p1[ 10:12, ] <- data.frame(extractPlotInfo(p1_safety), 3, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p1[ 13:15, ] <- data.frame(extractPlotInfo(p1_services), 2, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p1[ 16:18, ] <- data.frame(extractPlotInfo(p1_culture), 1, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p1[ 1:3, ] <- data.frame(extractPlotInfo(p1_pca_all), 6, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)

plot_mat_p1B <- data.frame(est = rep(0, 3), lower95 = 0, upper95 = 0, lower90 = 0, upper90 = 0, x = 0, D = NA, stringsAsFactors = F)
plot_mat_p1B[ 1:3, ] <- data.frame(extractPlotInfo(game_mod, 3:5), 1, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)

plot_mat_p2 <- data.frame(est = rep(0, 9), lower95 = 0, upper95 = 0, lower90 = 0, upper90 = 0, x = 0, D = NA, stringsAsFactors = F)
plot_mat_p2[ 1:3, ] <- data.frame(extractPlotInfo(p2_pca), 3, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p2[ 4:6, ] <- data.frame(extractPlotInfo(p2_neighbors), 2, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)
plot_mat_p2[ 7:9, ] <- data.frame(extractPlotInfo(p2_culture), 1, c("Common Humanity", "Countering Stereotypes", "Norms"), stringsAsFactors = F)


# Plotting ----------------------------------------------------------------

panelA <- ggplot(data = plot_mat_p1, aes(x = est, y = x, color = D, shape = D)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggstance::geom_linerangeh(aes(xmin = lower95, xmax = upper95), position = ggstance::position_dodgev(.6)) +
  ggstance::geom_linerangeh(aes(xmin = lower90, xmax = upper90), position = ggstance::position_dodgev(.6), size = 1.2) +
  geom_point(position = ggstance::position_dodgev(.6), size = 3, fill = "white") +
  xlab("Estimated Treatment Effect") +
  ylab(NULL) +
  scale_y_continuous(
    breaks = 6:1,
    labels = c(
      "Immigration Index\n(SD=1)",
      "Welcome Immigrant\nNeighbors (1-7)",
      "Support Increasing\nImmigration (1-7)",
      "Immigrants Increase\nNeighborhood Safety\n(1-7)",
      'Immigrants Contribute\nto Govt. Services\n(1-7)',
      "Immigrants Enrich\nCulture (1-7)"
    )
  ) +
  scale_x_continuous(breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, .8)) +
  coord_cartesian(xlim = c(-0.5, 0.85), ylim = c(.5, 6.5), expand = F) +
  scale_color_grey(name = "Treatment:", guide = "none") +
  scale_shape_manual(values = c(15, 21,17), name = "Treatment:", guide = "none") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), legend.position = "bottom") +
  ggtitle("A. Attitudinal Effects, Pilot 1 (N=351)")

panelC <- ggplot(data = plot_mat_p1B, aes(x = est, y = x, color = D, shape = D)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggstance::geom_linerangeh(aes(xmin = lower95, xmax = upper95), position = ggstance::position_dodgev(.5)) +
  ggstance::geom_linerangeh(aes(xmin = lower90, xmax = upper90), position = ggstance::position_dodgev(.5), size = 1.2) +
  geom_point(position = ggstance::position_dodgev(.5), size = 3, fill = "white") +
  xlab("Estimated Treatment Effect") +
  ylab(NULL) +
  scale_y_continuous(
    breaks = 1,
    labels = c("Difference in Tokens\n Given to Foriegn\n and Native Partners\n(0-10)")
  ) +
  scale_x_continuous(breaks = c(-0.4, -0.2, 0, 0.2, 0.4, 0.6, .8)) +
  coord_cartesian(xlim = c(-0.5, 0.85), ylim = c(.5, 1.5), expand = F) +
  scale_color_grey(name = "Treatment:", guide = guide_legend(reverse=TRUE)) +
  scale_shape_manual(values = c(15, 21,17), name = "Treatment:", guide = guide_legend(reverse=TRUE)) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), legend.position = "bottom") +
  ggtitle("C. Behavioral Effects, Pilot 1 (N=1,448)")

panelB <- ggplot(data = plot_mat_p2, aes(x = est, y = x, color = D, shape = D)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggstance::geom_linerangeh(aes(xmin = lower95, xmax = upper95), position = ggstance::position_dodgev(.25)) +
  ggstance::geom_linerangeh(aes(xmin = lower90, xmax = upper90), position = ggstance::position_dodgev(.25), size = 1.2) +
  geom_point(position = ggstance::position_dodgev(.25), size = 3, fill = "white") +
  xlab("Estimated Treatment Effect") +
  ylab(NULL) +
  scale_y_continuous(
    breaks = 3:1,
    labels = c(
      "Immigration Index\n(SD=1)",
      "Welcome Immigrant\nNeighbors (1-7)",
      "Immigrants Enrich\nCulture (1-7)"
    )
  ) +
  scale_x_continuous(breaks = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, .8)) +
  coord_cartesian(xlim = c(-0.5, 0.85), ylim = c(.5, 3.5), expand = F) +
  scale_color_grey(name = "Treatment:", guide = guide_legend(reverse=TRUE)) +
  scale_shape_manual(values = c(15, 21,17), name = "Treatment:", guide = guide_legend(reverse=TRUE)) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), legend.position = "none") +
  ggtitle("B. Attitudinal Effects, Pilot 2 (N=737)")


# Panel combination and saving --------------------------------------------

Dlegend <- get_legend(panelC + guides(color = guide_legend(title.position = "top", reverse = T), shape = guide_legend(title.position = "top", reverse = T)) + theme(legend.title = element_text(hjust = 0.5), legend.margin = margin(30, 0, 0, 0)))
plots <- align_plots(panelA, panelC + theme(legend.position = "none"), panelB, Dlegend, align = 'vh', axis = "lt")
left <- plot_grid(plots[[1]], plots[[2]], ncol = 1, rel_heights = c(1, 2/6))
right <- plot_grid(plots[[3]], plots[[4]], ncol = 1, rel_heights = c(1, 2/6))
fig_s8 <- plot_grid(left, right, ncol = 2)

ggsave(here("SI", "figure_s8.pdf"), fig_s8, height = 8.5, width = 11)

